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ABSTRACT 

^—i We analyse three-dimensional smoothed particle hydrodynamics simulations of the 

Oh evolution of initially quasi-circular massive black hole binaries (BHBs) residing in the 

central hollow (cavity) of massive self-gravitating circumbinary discs. We perform a 
set of simulations adopting different thermodynamics for the gas within the cavity and 
for the 'numerical size' of the black holes. Aimed at understanding the physical roots 
of the secular evolution of the binary, we study the interplay between gas accretion and 
gravity torques in changing the binary elements (semi-major axis and eccentricity) and 
its total angular momentum budget. We pay special attention to the gravity torques, 
by analysing their physical origin and location. We show that (i) the BHB eccentricity 
grows due to gravity torques from the inner edge of the disc, independently of the 
accretion and the adopted thermodynamics; (ii) the semi-major axis decay depends 
not only on the gravity torques but also on their subtle interplay with the disc-binary 
angular momentum transfer due to accretion; (iii) the spectral structure of the gravity 
torques is predominately caused by disc edge overdcnsitics and spiral arms developing 
in the body of the disc; (iv) the net gravity torque changes sign across the BHB 
corotation radius: gas inside this radius exerts a net positive torque, while streams 
located outside this radius (but within the cavity) exert a net negative torque. The 
relative importance of the two might depend on the thermodynamical properties of 
the instrcaming gas and is crucial in assessing the disc-binary angular momentum 
transfer; (v) the net torque manifests as a purely kinematic effect as it stems from 
the low density cavity, where the material flows in and out in highly eccentric orbits. 
Thus both accretion onto the black holes and the interaction with gas streams inside 
the cavity must be taken into account to assess the fate of the binary. 

Key words: Black hole physics - Accretion, accretion discs - Numerical - Hydro- 
dynamics 



1 INTRODUCTION 

In the curre ntly favoured h i erarch ical framework of structure 
formation l)White fc Reed Il978t ). galaxies evolve through 
a complex sequence of merger and accretion events, and 
the existence of massive black holes (BHs) at their cen- 
tres is nowadays a w ell established observational fact (see 
iGiiltekin et alj [20091 . and references therein). By combin- 
ing together these two pieces of information, the forma- 
tion of a large number of massive black hole binaries 



(BHBs), following galaxy mergers throughout cosmic his- 
tory, is a natural consequence of the structure formation 
process (|Begelman et alJll98ch . Although this is corrobo- 
rated by several o bserved quasar pairs at a ~ lOOkpc pro- 
jecte d separation (iHennawi et al.l I200C ; iMvers et alj 120071 . 



120081 ; iForeman et al.l 120091 ; IShen et alJboilf ). and by few 
< k pc dual accreting BHs embedded in the s ame galaxy 
(e.g. iKomossa et al1l2003l ; iFabbiano et al.ll201lh . identifica- 
tion of gravitationally bound BHBs in galaxy centres re- 
ma ins elusive (for an up-to-date review on the candidates 
see lDotti et alj|2012l and references therein). 
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their existence, much theoretical work has focused lately on 
Keplerian massive BHBs residing in galactic nuclei. One rea- 
son is that a deep understanding of the interplay between 
BHBs and their dense (stellar and gaseous) environment is 
required to predict robust signatures that may allow their 
identification. Additionally, it is still unclear how the gap be- 
tween the two theoretically well understood stages of BHB 
evolution is bridged: (i) the dynamical friction driven stage, 
when the two BHs spiral in toward the centre of the merger 
remnant down to pc separations and (ii) the final inspiral 
driven by gravitational waves (GWs), which become effi- 
cient when the two BHs are at a separation <J10 -2 pc. Both 
dense stellar and gaseous environments have been shown 
to be effective in extr acting the binary energy and angular 
momentum (see, e. g . , lEscala et ai1l2005l; iDotti et al.l 120071 : 
ICuadra et all 120091 ; iKhan et all l201ll ; IPreto et al.1 boilf ). 
likely driving the system to final coalescence (an extensive 
discussion on the fate of sub-parsec BHBs can be found in 
IDotti et al.1 120121 ). Scenarios involving cold gas are partic- 
ular appealing not only because they might produce dis- 
tinctive observational signatures, but also because cold gas 
dominates the baryonic content in most galaxies at redshifts 
higher than one, providing a natural reservoir of energy and 
angular momentum to drive the BHB towards coalescence. 

Numerical simulations of wet galaxy mergers (e.g., 
iMihos fc Hernauisd[l996l ; lMaver et alj|2007h have led to the 
following picture for the post merger evolution. Cold gas 
is funnelled to the central w 100 pc by gravitational in- 
stabilities, where it forms a puffed, rotationally supported 
circumnuclear disc. Disc-like structures are actually ob- 
served in ULIRGs which are tho ught to be gas rich post- 



merger star-forming galaxies (e.g., Sanders fc Mirabel 19961; 
Downes fc Solomorll998l ; lDavies et al.ll2004al lbl; lGreve et al l 



2009). In the models, the two nuclear BHs efficiently spi- 



ral to sub-pc scales owing to dynamical friction against 
the massive circumnuclear disc (|Dotti et al.l l2007h . eventu- 
ally opening a cavity (or h ollow) in the gas distribution 
l|Goldreich fc Tremaindll980l) . The subsequent evolution of 
the system is determined by the efficiency of energy and an- 
gular momentum transfer between the BHB and its outer 
circumbinary disc. 

The investigation of coupled disc-binary systems has a 
long-standing tradition in the context of planetary dynamics 
llGoldreich fc Tremame1ll980l;lLin fc Papaloizoulll986l; IWardl 
19971 ; iBrvden et al.lll999l ; iLubow et al.lll999l ; iNelson et al.1 
2000), where the focus usually lies on the extreme mass ratio 
situation, i.e., a star surrounded by a circumstellar disc, with 
a planetary companion embedded in it. The comparable 
mass case has also been extens ively investigated in the con- 
text of binary star formation ijArtvmowicz fc Lubowlll994 
Bate fc Bonnei]|l997l;lGiinther fc Klevll2002l; iGiinther et al l 



2004), where a boost of activity has been triggered by 



imaging of nearby you ng binary stars em bedded in hol- 
low circumbinary discs l|Dutrev et al.lll994T ). More recently, 
the techniques adopted in these fields have been applied 
to the BHB case. In the last decade, several investi- 
gations were devoted to the study of comparable mass 
BHBs evolving in circumbinary discs, explo iting a vari- 
ety o f analytical and numer i cal te c hniqu e s lllvanov et al. | 
1999; Armitagc & Nataraiar J |2002l. |2005| ; Hayasaki ct al. 
2007 1. 120081; iMacFadven fc Milosavlievidl 120081; lHavasakJ 
20091 ; ICuadra et al l 120091 ; lLodato et all 120091 ; iNixon et all 



l201ll ; iRoedig et all l201ll ; IShi et al.1 l201ll ). However only 
few of them focu sed on the details of the dynamic al disc- 
binary interplay. MacF adven fc Milosavlievid (|2008T ) (here- 
inafter MM08) made use of two-dimensional grid-based hy- 
drodynamic al simulation to study a B HB embedded in a 
thin a-disk (|Shakura fc Sunvae v Hl973T ). They showed that 
the gas flowing through the cavity increases the energy and 
angul ar momentum t ransfer between the disc and the bi- 
nary. [ShTet^al] f|201 lh (hereinafter Sll) confirmed that re- 
sult using full 3D magnetohydrodynamics (MHD) simula- 
tions. However, in both studies the BHB was on a fixed cir- 
cular orbit, the central region of the disc (i.e. within twice 
the binary separation) was excised from the computational 
domain, and the disc self-gravity was neglected. 

We employ full 3D smoothed particle hydrodynamical 
(SPH) simulations to study the interaction between BHBs 
and their surrounding discs. Our goal is to give a detailed 
description of the coupled disc-binary dynamics, paying par- 
ticular attention to the competing effects of disc-binary 
gravitational torques and gas accretion onto the BHs in the 
evolution of the binary angular momentum budget. We si- 
multaneously evolve the disc and the two BHs in a self- 
consistent way. This enables us to separately investigate the 
effect of gravitational torques coming from different disc 
regions and to directly link different physical mechanisms 
(gravitational torques and accretion) to the evolution of in- 
dividual quantities describing the binary (mass, mass ra- 
tio, semi- major axis and eccentricity). Our approach allows 
a broader understanding of the coupled dynamics with re- 
spect to previous studies (MM08, Sll), where the binary 
was modelled as a fixed forcing quadrupolar potential and 
the central region was excised. Thanks to the relative low 
computational cost of the SPH implementation we could 
perform different runs varying physical and numerical pre- 
scriptions that may play a relevant role in the disc-binary 
energy and angular momentum exchanges. Specifically, we 
checked the possible dependence of our results on the 'nu- 
merical size' of the two BHs (i.e., their sink radii, see next 
Section), and on the adopted equation of state (EoS) for the 
gas within the central cavity. 

The paper is structured as follows: we first describe the 
numerical setup and initial conditions in Section [2J then 
we show the importance of both accretion and gravitational 
torques on the binary evolution in Section[3] We study in de- 
tail the origin and strength of the mutual disc-binary gravi- 
tational torques in Section|3]highlighting the influence of the 
thermodynamics prescription. In Section[5]we briefly discuss 
the temporal behaviour and the spatial geometry of the ac- 
cretion, speculating on the evolution of the binary mass ratio 
and BH individual spins. Finally, we compare our results to 
previous work in Section HJ] and discuss the implications of 
our findings and give our conclusions in Section [7] 



2 SIMULATIONS 

The model and numerical setup of this work is closely relate d 



to that of IRoedig et~al] i|201ll ) and ICuadra et all (2009), 



hence we only outline the key aspects in the following, and 
refer the reader to these two papers for further details. 

We simulate a self-gravitating gaseous disc of mass 
Md = 0.2 M around two BHs of combined mass M = 
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Figure 1. Relaxed meridional density maps of the disc for the 
adia (top) and iso (bottom) runs; the black dots indicate the 
BHs, the axis are in units of the binary semi-major axis a. Note 
that this is not an azimuthal average but a vertical slice of the 
disc. 
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Figure 2. Time and azimuth-averaged values of the disc scale- 
height H/r and the disc eccentricity e^isc as a function of radius 
for the different simulations. 



Mi + M%, mass ratio q = M2/M1 = 1/3, eccentric ity e and 
semi- major axis a, using the SPH-code Gadget-2 (jSpringell 
120051 ) in a modified version that in cludes sink-particles 
which model accret ion on to the BHs fSp ringel et al.l [20051 ; 
ICuadra et ai1l2006t ). Moreover, the orbit of the BHB is fol- 
lowed very accurately by using a fixed small time-step and 
summing up directly the gra vitational force from every other 
particle in the simulation l|Cuadra et al.l 12009] ). The disc, 
which is co-rotating with the BHs, radially extends to about 
7a, contains a circumbinary cavity of radial size ~ 2a and is 
numerically resolved by 2 million particles. The numerical 
size of each BH is denoted as r s i n k, the radius below which 
a particle is accreted, removed fr om the simulatio n, and its 
momentum is added to the BH l|Bate et al.lll995l ). The gas 
in the disc is allowed to cool on a time scale proportional to 



Table 1. Initial parameters and names of the runs in unitlcngths 
of the code. 



Model 


^sink 


Cavity EoS 


Disc EoS 


adialO 


0.1 


adiabatic+/3 


adiabatic +/3 


isolO 


0.1 


isothermal 


adiabatic +/3 


adia05 


0.05 


adiabatic+/3 


adiabatic +/3 


iso05 


0.05 


isothermal 


adiabatic +/? 



the local dynamical time of the disc tdyn = fo 1 — 27r/f2o, 
where fio = (GMo/oo) 1 ' 2 is the initial orbital frequency of 
the binarjQ. To prevent it from fragmenting, we force the 
gas to cool slowly, se tting f} = t coo \/t AYn = 10 l|Gammiel 
l200ll ; lRice et alj|2005h . We use two different treatments for 
the thermodynamics inside the cavity, denoted as adia and 
iso, respectively. In the adia runs, the gas within the cav- 
ity is allowed to both cool via the /3 prescription and to 
heat up adiabatica lly, as it is the remainder of the disc 
l|Cuadra et al.| [2009). In the iso runs, we define a threshold 
radius r cav i ty = 1.75a, below which the gas is treated isother- 
mally, meaning that the internal energy per unit mass is set 
to be u « 0.U(GM/R) l|Roedig et al.ll201ll ). 

As initia l condi tions we use a relaxed snapshot from 
Cuadra et al.l 



taken at their t 



500fin 



(see 
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Section 2.1J0 . The nomenclature of 
the runs and the parameters used are listed in Tab.[T] Each 
run is evolved for » 90 binary orbits, and we store the out- 
put in single precision ~ 6 times per orbit. We also closely 
track accretion by storing the position and velocity of the 
two BHs and of each accreted particle at the time the latter 
crosses the sink radius. For both prescriptions (adia and iso) 
we show the relaxed density configurations sliced along the 
meridional plane in Fig. Q] (see also Figs. [8}{9] for a face-on 
viewjfl and the measured, time-averaged values for both the 
scale height of the disc H/r, and the eccentricity of the disc 
edisc in Fig. [2] H is taken to be the height above and below 
the disc-orbital plane in which 70% of the mass inside the 
annulus (r, r + dr) is found. In their bodies (i.e. at r > 2a), 
the discs may be considered thin, as shown in the upper 
panel of Fig. [2] with a maximum thickness H/r « 0.2, and 
quasi circular (lower panel of Fig. [2j . Generally, we can con- 
sider the physics to be resolved if the smoothing length h is 
smaller than the characteristic lengthscale: radially the cri- 
terion h/r <C 1 is achieved well: h/r G (0.005,0.2); whereas 
vertically the resolution in the main body of the disc is about 
h/H G (0.1,0.25)0. 

Unless otherwise stated, we will present all the relevant 
quantities in the natural units of the simulation by setting 
G — Mo = ao = 1. It follows that also the initial circular 
velocity of the binary is V c ,o = 1 and its initial period Po = 
2n. We will then discuss the astrophysical implications of 



1 Subscripts refer to the initial values of any parameter. 

2 Notice that in the remainder of this article we refer to the time 
of this snapshot as t = 0. 

3 Plots made using SPLASH jPrice!l2007ll . 

4 The specified values are not met inside the cavity, due to the 
relatively small quantity of gas. This yields h/H G (1,4) for r < 
1.5a. 
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Figure 3. Semimajor axis (top panel in each plot) and eccen- 
tricity (bottom panel in each plot) evolution for all the four runs. 
In each panel, the red line represents the full evolution directly 
taken from the simulation whereas the black line is the equivalent 
evolution when considering the energy and angular momentum 
exchanges due to accreted particles only. 



our findings by scaling them to fiducial astrophysical BHB 
systems. 



3 BINARY EVOLUTION: CAUSES 

For each of the four runs, we measure the two primary or- 
bital elements: eccentricity e and semi-major axis a. Their 
evolution is shown by the red lines in Fig. [3] In all four runs, 
we find a(t) < and e(t) > 0. While the eccentricity evo- 
lution is largely independent on the sink radius value and 
the adopted EoS within the cavity, the orbital shrinking is 
much faster in the adia runs, in which the binary shrinks 
by 4 — 5% over 90 orbits. Conversely, in the iso runs, the 
two BHs get only about 1% closer. Linear extrapolation of 
such results at face value implies binary coalescence on a 
timescale of ~ 2000Po and ~ 10 4 Po for the adia and the iso 
runs, respectively. Whereas, assuming a const ant eccentric- 



ity_gj 

|20T 



rowth rate, the limiting value predicted bv lRoedig et al.l 



2011) would be reached in less than ~ IOOOPq. 



3.1 Gravitational torque and accretion 

contribution to the angular momentum 
budget 

Being interested in the physical mechanisms driving the bi- 
nary evolution, we firstly identify two distinct processes: 

(i) the gravitational torques exerted by the gaseous par- 
ticles onto each individual BH, Tg (gravity torque) 

(ii) the accretion of instreaming particles crossing either 
BH sink radius, (dL/di) acc (accretion torque). 



Conservation of the total angular momentum in the simula- 
tions implies 



dL d~L 

at dt acc 



(1) 



where L is the BHB orbital angular momentum vector 0. All 
vectors are computed with respect to a Cartesian reference 
frame centred in the BHB centre of mass (CoMfl the bi- 
nary initially lies in the x-y plane with angular momentum 
oriented along the positive z axis. In our SPH simulations, 
Tg can be computed at each snapshot by direct summation 
of individual particle torques onto each BH yielding 



T G 



N 2 



GM k m,j(rj - r k ) 
|r, - r fc | 3 



(2) 



where j runs over all N gas particles, k identifies the two 
BHs, r are position vectors, and m and M denote particle 
and BH masses respectively. On the other hand, (dL/di) acc 
can be computed by assuming instantaneous linear momen- 
tum conservation of each accreted particle yielding AL = 
rfe x rrijVj. Here r& is the position vector of the accreting 
BH at the moment of swallowing the particle, and v,- is the 
particle velocity vector. Numerically, we can evaluate the 
binary angular momentum change AL at each snapshot by 
writing Eq. fT} in the form 



AL = T G At + y~] r k x m 3 v 3 - 

At 



(3) 



where At is the interval between two subsequent snapshots, 
and the sum runs over all the accretion events occurring in 
this time lapse. 

We use Eq. © to assess the relative importance of grav- 
itational torques and accretion in the evolution of the sys- 
tem. In Fig. 3] we plot the evolution of the x, y, z components 
of the angular momentum L as directly evaluated from the 
positions and velocities of the two BHs stored in the simula- 
tion outputs (red), together with the relative changes due to 
accretion (black dotted) and gravity torques (black dashed) 
according to Eq. ([3)l . The solid black line is the sum of the 
latter two which should overlap with the red line if our de- 
composition is sufficiently accurate and these are the only 
two physical processes determining the binary evolution. For 
comparison, in simulation units, the initial binary angular 
momentum is ~ 0.18. The L x and L y components are al- 
most constant, showing fluctuations at the 10~ 4 level. The 
mismatch between the red and the solid black lines here is 
because simulation outputs are in single precision, i.e. ac- 
curate to the fourth digit, which is the fluctuation level of 
the computed quantity. On the contrary, the L z component 
changes up to a 10~ 2 absolute level (i.e., about 5% of the 
initial value), and in this case we see very good agreement 
between the two lines. Note that in three of the runs, the 
overall angular momentum grows over time, even though a 
shrinks and e increases. This counter intuitive result is due 



5 Throughout the paper we will denote all vectorial quantities in 
boldface. 

6 The binary CoM slightly wiggles around the total binary-disc 
CoM. As a cross check, we also computed all the relevant quanti- 
ties with respect to the latter, finding no significant difference in 
the results. 
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Figure 4. Consistency check of the evolution of the binary angular momentum components for all the simulations. In each plot, the L x , 
L y and L z component evolution is shown from the top to the bottom. In each panel, the dashed line is the AL exchange balancing the 
gravity torques at each time step, the dotted line is the AL exchange computed from the accreted particles, and the solid black line is 
the sum of those two. For comparison, the red line is the angular momentum evolution taken directly from the raw simulation data. All 
AL are in units of [MoaoV c ,o]- 



to the dominance of accretion in the angular momentum 
budget, and will be explained in § 13.21 below. 

As a general rule, the accretion contribution to the bi- 
nary L evolution (i.e., the accretion torque) is at least com- 
parable to the effect of the gravitational torques. It is there- 
fore also interesting to see the accretion contribution to the 
evolution of the binary orbital elements. Having stored all 
the accretion events, this can be done separately by simply 
evolving the binary from ao and eo only by adding one ac- 
creted particle after the other, imposing linear momentum 
conservation. This "would-be" evolution of e(t) and a(t), ac- 
counting only for the accretion effect, is shown in black in 



all panels of Fig. [3] It is clear that e(t) is mainly driven 
by gravitational torques, with accretion playi ng a negligible 
role. This is in line with our previous findings l|Roedig et al.l 
l201ll ). where the evolution of the eccentricity was attributed 
mainly to the gravitational interaction between the binary 
and the overdensities excited by its quadrupolar potential in 
the inner rim of the circumbinary disc. Conversely, the a(t) 
plots highlight the importance of accretion for the semimajor 
axis evolution. In the iso05 case, the binary shrinking due to 
accretion is actually larger than the total one, meaning that 
the disc torques would force t he binary to expand; an ef- 
fect similar to that described bv lLin fe Papaloizoul (|201ll ) in 
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the context of planetary migration in self-gravitating discs. 
Note that in the adia cases a is dominated by the effect of 
the disc. This explain s the m atch between such simulations 
and the llvanov et al.l (|l999r ) analytical model based on an- 
gular momentum transport through the disc (|Cuadra et al.l 
l2009h . 



3.2 Dissection of the binary evolution into its 
components 

So far, we have separated the relative contribution of grav- 
itational torques and accretion to the binary angular mo- 
mentum budget. We now investigate how the angular mo- 
mentum change is distributed among the relevant binary 
quantities. As clearly shown in Fig. 01 L x and L y show only 
small fluctuations (and therefore remain small compared to 
L z , since the binary initially orbits in the x-y plane). From 
here on, we will therefore concentrate on the dominant L z 
component. The binary angular momentum is 



L z = ^GMa(l-e 2 ), 



(4) 



where fi — M1M2/M is the binary reduced mass. Eq. 
can be differentiated to separate the contribution of each 
single relevant quantity to the angular momentum change: 

L z a M (i e . . 

T z = Ya + 2M + £ ~ T~^ e ' ( ) 

By directly measuring a, , M , from the simulation, we can 
evaluate each single term and decompose the L z evolution 
accordingly. This is shown in Fig. [5] for all the four runs. 
Note that the sum of the four contributions (solid black line 
in each panel) closely matches the measured L z evolution 
directly measured from the simulation (red), validating our 
first-order expansion. 

It is worth noting that the increase of L z is perfectly 
compatible with semi-major axis shrinkage and eccentricity 
growth. In fact, Eq. ([5)l can be inverted to obtain 



a _2T Z M 2fi 

a ~ T7 ~ M ~ + I 



2e 



(6) 



where we used the fact that L z = T z . Even if the total 
torque is positive, and the eccentricity grows, the binary 
can shrink because of the negative contribution of the M 
and the /x terms. In particular, as we will see below, the 
higher accretion onto the secondary hole results in a large 
increase of n (long-dashed line in Fig. [5]) which is the main 
driver of the binary angular momentum growth. Eqs. ([5]) 
and ([6]) highlight the importance of accretion in determining 
the binary temporal evolution. 



4 GRAVITATIONAL TORQUES 

As pointed out in the previous Section, understanding sec- 
ular dynamics of BHBs in gaseous environments is closely 
tied to studying the interplay between long distance gravita- 
tional forces, hydrodynamics and accretion. In this Section, 
we focus on the torques coming from the self-gravitating 
gaseous disc, investigating the influence of regions at differ- 
ent distances from the BHB. 
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Figure 5. Decomposition of the L z evolution for the four runs. 
The AL Z contributions stem from a (dotted), e (dashed), M (dot 
dashed) and fj, (long dashed) variations separately. The sum of 
the four contributions is given by the solid black line, whereas the 
solid red line is, as in Fig. [4] the angular momentum evolution 
taken directly from the raw simulation data. All AL are in units 
of [M a V c ,o]. 



4.1 Time averaged torque profiles 

Gravitational torques exerted by the disc onto the BHB can 
be directly calculated at each snapshot from Eq. ((2j) . It is of 
great interest to understand where such torques originate, 
and given the cylindrical nature of the problem, it is natural 
to investigate their azimuthally-averaged radial distribution. 
We take r to be the projected radial coordinate (in the x-y 
plane) from the binary CoM, and we define dT/dr to be the 
differential torque, integrated over the azimuthal coordinate 
and the disc height. The total average torque exerted by gas 
in the projected distance interval [a, b] from the binary CoM 
is 



(7kb,) = (/ ^ 



(7) 



where (•) denotes temporal averaging over the entire simu- 
lation. 

In Fig. [6] we show the time averaged differential torque 
acting on the binary, (dT/dr) (blue), together with its in- 
tegral according to Eq. Q (black). We also decompose the 
former in the two components acting on the primary (green) 
and the secondary (red) BH. All torques are null at the bi- 
nary corotation radius we therefore show the total torque by 
integrating from this point inwards ((Tm qi), binary region) 
and outwards ((Tj li00 ]) ; disc region). 

In all simulations, the local average torque shows an 
oscillatory behaviour with a sharp maximum at the loca- 
tion of the secondary BH, r ~ 0.75, and a deep minimum 
in the cavity region at 1 < r < 2. In the body of the disc 
(r > 2) positive and negative peaks alternate, but they are 
shallow and almost cancel out, giving a negligible contribu- 
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Figure 6. Differential torques dT/dr in units of [GMq<2<7 2 ] aver- 
aged over the entire simulations. In each panel we show the differ- 
ential torque on the primary (green), on the secondary (red), the 
sum of the two (blue), and the total integrated torque (black) ac- 
cording to Eq. [7] This latter is integrated starting from a inwards 
and outwards. Notably, the inward torque is positive, whereas the 
outward torque is negative. 



tion to the total torque (as witnessed by the fact that the 
integrated torque is basically constant for r > 2). Note that 
the torque on the secondary BH is always larger than the 
torque on the primary, due to its proximity to the outer disc 
resulting in a stronger interaction. The general behaviour is 
qualitatively the same in the adia and iso runs. In this lat- 
ter case, however, we found a much sharper negative peak 
at r ~ 1.7 followed by a smaller secondary negative bump 
at r ~ 1.2. This appears to be an artefact of the sudden 
change in EoS of the gas inside the cavity, at r = 1.75. The 
net result is a smaller negative (Tri iOC ]) which has important 
consequences on the binary evolution. We see in fact that in 
the iso runs (Tr li0 ]) and (Tn i00 i) almost cancel out, meaning 
that, overall, gravitational torques do not change the binary 
angular momentum. Conversely, in the adia runs (Tn l) is 
much smaller (in absolute value) than (Tn i00 ]), implying an 
efficient ang ular momentum tra nsfer from the binary to the 
gas (MM08; ICuadra et alj|200sh . 



4.2 Spectral analysis 

We now consider the torque evolution in time. We separately 
discuss torques coming from different disc regions by show- 
ing both the time series and their associated power spectra. 
We cut the spatial domain into the following radial annuli: 



(i) 

(ii) 

(iii) 

(iv) 

W 



a < r < 10 a 

a < r < 1 a 

1 a < r < 1.8a 
1.8a < r < 2.5a 
2.5a < r < 10 a 



the entire domain 
the 'binary region' 
the 'cavity region' 
the 'rim region' 
the 'disc region' 



The associated time series and power spectra are shown 
respectively in the left and right panels of Fig. [7] 

In all the simulations, the overall torque (i), shows a 
clear periodic oscillation, much larger in amplitude than its 
average value. The power spectrum unveils several distinc- 
tive peaks, with relative amplitudes that can vary signifi- 
cantly for different simulations. In particular, peaks in the 
iso simulations are much sharper and better defined. This 
is because in the adia simulations the BHB shrinks signif- 
icantly, resulting in a broadening of the characteristic fre- 
quencies. Moreover, as we shall see in the next section, the 
disc sub-structures are much better defined in the iso runs, 
giving rise to neater features. 

In the binary region (ii) the torque is mostly coherent 
and positive. Because the torque strength in (ii) is regulated 
by the mass inflow, it shows periodicities that are related to 
the accretion flow: a disc component around / ~ 0.25P _1 
(corresponding to the disc peak density at r ~ 2.5a), the 
forcing frequency of the binary at / = P^ 1 , and the beat 
between thes e two at / w 0.75P _1 (see Section [5] and 
iRoedig et al.l (|201ll ) § 5.1). In the cavity region (iii) the 
torque is negative on average but strongly oscillating. Sev- 
eral periodicities are detectable, the most striking being a 
peak at / w P _1 which appears again to be directly related 
to the binary period. In the rim region (iv) the torque is 
highly oscillating, and the strongest feature is a sharp peak 
at / ~ 1.3P _1 , whereas in the disc region (v) the only sig- 
nificant spectral component is at / m 1.7P^ ■ As a general 
trend, moving from the inner region to the disc body, torques 
become incoherent (i.e. they average to zero) and strongly 
oscillating (compare the power spectra scales in the different 
panels of each plot). 



4.3 Interpretation: torque origin and location in 
the disc 

In this Section, we provide a global interpretation of the 
features observed both in the radial distribution of the time 
averaged torques and in their temporal evolution. The ar- 
guments discussed below are supported by Figs. [8n9l 1101 and 
by Tab. d 



4-3.1 Origin of the positive and negative peaks 

It seems natural to compare the torque radial profiles ob- 
tained by our simulations with linear theory perturbation, 
in which torque minima and maxima are connected to outer 
Lindblad resonances (OLRs). It is in fact tempting to asso- 
ciate the torque minimum at r fn 1.6a with the 2:1 OLR. 
We should however be careful in pushing this interpretation 
too far, since, as already pointed out by MM08 and Sll, 
the assumptions of linear theory are not satisfied in this 
context. Most importantly, looking at the upper panels of 
Figs. [8] and [9] we notice that the region r < 2a is almost de- 
void of gas, and the streams are almost radial. Mass fluxes 
reported in Tab. clearly show that, at any radius, there 
are always fluxes of ingoing and outgoing mass, resulting in 
a steady net inward flux consistent with the accretion rate 
onto the two BHs. A strict OLR interpretation of the torque 
minimum at r ~ 1.6a would instead require particles in cir- 
cular orbit at that radius, experiencing a secular effect due 
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Figure 7. Each plot depicts the torques exerted by the disc on the BHB in the time domain (left panels) and in the frequency domain 
(right panels). In each plot, from the top to the bottom, pairs of panels refer to the five domains discussed in the text: total torque (i), 
torque exerted by the material located at r < a (ii), a < r < 1.8a (iii), 1.8a < r < 2.5a (iv) and r > 2.5a (v). In each left panel, the 
red line is the raw oscillating torque, and the blue one is the torque smoothed over three periods to show the average behaviour. In the 
associated right panel, we plot the power spectrum as a function of frequency in units of the initial binary orbital frequency. All torques 
strongly oscillate, showing several characteristic frequencies, however, note the different scales of the panels. 



to the phase-coherent periodic forcing of the binary; this is 
not what happens within the cavity region. The strongest 
2:1 OLR is certainly responsible for the evacuation of the 
gas close to the binary and for the formation and mainte- 
nance of a cavity, however cannot be directly responsible for 
the coherent torque seen in the cavity region. This is also 
supported by the fact that MM08 and Sll find a minimum 
at the same location (r « 1.6 — 1.7a) for an equal mass 
binary, where the 2:1 resonance is absent (because of the 
symmetry of the forcing potential), and the location of the 
strongest OLR (3:2) would be at r w 1.3a. The strong neg- 



ative torque in the cavity region has a purely kinetic origin: 
material ripped off the disc edge forms well defined streams 
following the two BHs, which are clearly distinguishable in 
both the surface density plots shown in Fig. [S] and [5] The 
streams are responsible for the yellow tails following the two 
BHs at ~ 1.5a in the torque density panels, which lead to 
a net negative torque. Conversely, at r,>2a, we have a well 
defined, almost circular disc, and the torque density peaks 
at r « 2a and r w 2.5a ca n be identified with the loci o f the 
strong 3:1 and 4:1 OLRs l|Artvmowicz fc Lubow|[l994l ). 

Our simulations also allow us to investigate the torques 
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Figure 8. Top panel: typical instantaneous logarithmic surface 
mass density of the iso05 simulation. Bottom row: combined sur- 
face torque density exerted by the disc onto the binary (left 
panel), onto the primary BH (central panel) and onto the sec- 
ondary BH (right panel). All plots are in code units: G = Mo = 
a = 1. 




Figure 9. Same as Fig. [8] but for the adia05 run. Here the fea- 
tures highlighted above are even more evident (especially the sym- 
metric overdensities at the inner edge of the disc). 
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Figure 10. Disc properties: average surface density (top panel), 
and Toomre parameter Q (bottom panel) as a function of r. 
Line style denotes isolO (long dashed), iso05 (dot dashed), adialO 
(solid) and adia05 (dashed). 



within the binary corotation radius at r < a, a region of- 
ten excised in grid-based simulations (see MM08 and Sll). 
Here we find strong positive torques on both BHs. This is 
because the infalling material approaches the BHs at super- 
Keplerian velocities, and bends in a horseshoe fashion, ex- 
erting a net positive torque in front of them. In fact, the 
maximum positive torque basically coincides with the loca- 
tion of the two BHs (sharp peak at r ~ 0.75a for the sec- 
ondary and a broader peak around r ~ 0.3a for the primary, 
see Fig. [6]). The very same effec t, in the context of plan e- 
tary migration, was discussed bv lLin fc Papaloizoul (|201ll )[1 
Note that the positive torque is related to this 'stream bend- 
ing', and not directly to the small discs of gas orbiting each 
BH; its magnitude is in fact similar in the iso and in the 
adia simulations, even though in the former, the mass in the 
minidiscs is significantly larger as shown by the time and 
azimuthally averaged surface density profiles in the upper 
panel of Fig. 1101 



4-3.2 Disc structure and characteristic torque frequencies 

The surface density profiles in Fig. [8] and [9] highlight some 
clear differences between the iso and the adia runs: (i) the 
appearance and shape of the minidiscs; (ii) the amount of 
gas feeding the streams; (iii) the definition of the inner disc 
edge. In the iso runs, we find a large amount of gas forming 
mini-discs around both BHs and an almost empty cavity 
except for two tenuous streams connecting the outer edge 



7 iLin fc Papaloizou (2011) found that such positive torque can 
in principle counter-balance the negative torque due to the disc 
and the tails, causing, in their case, angular momentum to be 
transferred to the planet, resulting in an outwards migration. 



10 C. Roedig et al. 



Table 2. Accretion rates onto the binary and mass fluxes F in units M/Pq x 10 -5 at two selected distances to the binary CoM: n = a 
and T2 = 1.5 a. Mi and M2 are the average accretion rates on M\ and M2 respectively, while M is the sum of the two. At both selected 
distances, Fj n and -F out represent the ingoing and outgoing fluxes, and F net = F; n — -F out is the net ingoing flux. All quantities have been 
averaged over 90 binary orbits. 













n = a 






T2 = 1.5a 




Model 


Mi 


M 2 


M 






Fnetx 




-Fout 2 


-Fnct 2 


isolO 


4.5 


10.7 


15.2 


19.6 


4.5 


15.1 


33.7 


18.8 


14.9 


iso05 


1.2 


9.2 


13.4 


19.1 


5.5 


13.6 


31.7 


18.4 


13.4 


adialO 


3.8 


9.5 


13.3 


17.6 


4.3 


13.3 


42.8 


29.7 


13.2 


adia05 


3.2 


1.6 


7.8 


13.5 


5.8 


7.7 


36.0 


28.4 


7.6 



to the mini-discs. The edge of the disc is quite circular and 
well defined. Conversely, due to the gas inside the cavity be- 
ing hotter in the adia runs, the streaming activity is more 
violent with thick spirals that do not settle into well de- 
fined mini-discs, but rather form a diluted, oo-shaped cloud 
around the binary. The edge of the disc is strongly disturbed 
by the ripped out streams and is usually not circular. The 
different streaming activity is also confirmed by numbers in 
Tab. [2] where we collect the average inwards and outwards 
mass fluxes at r = 1.5a and r = a. Mass fluxes are generally 
larger in the adia case. Note that this does not necessarily 
result in larger accretion, as outgoing fluxes are also larger 
in this case. Most of the material in the streams, does not 
end up in an accretion flow, but is accelerated back to the 
disc (a sort of 'slingshot') impacting on the inner edge, an 
effect already seen and extensively described by MM08 and 
Sll. The amount of impacting gas is 50% larger in the adia 
case, possibly contributing to the inner edge destabilization 
and to the larger perturbation in the disc structure. The 
appearance of spirals in the disc is related to the behaviour 
of the Toomre Q parameter, shown in the bottom panel of 
Fig. I lOt which quantifies the degree of self-gravity of a disc. 
Not surprisingly the disc develops a spiral pattern in the re- 
gion 3a <J r <J 5a, where, in fact, we find that t he Toomre pa- 
rame ter reaches its minimum value Q sa 1.5 (jCuadra et al.l 
2009). The shape of the spiral arms varies much in time, 
and their definition is different from simulation to simula- 
tion. However, at any time we find spiral structures with 2, 
3 or 4 arms, which remain confined between 3a < r < 5a. 

The disc structures highlighted above are reflected in 
the torques depicted in the lower panels of Fig. [S] and [5] The 
total torque shows a typical quadrupolar structure, which 
leads to rapid oscillations averaging to zero in the disc body. 
Conversely, as already discussed, in the inner cavity we can 
appreciate the yellow tails following the two BHs at~ 1.5a, 
leading to a net negative torque. 

The main peaks observed in the torque power spectra 
must therefore arise from the structures we just described. 
In particular we identified the overdensities at r ~ 2a in 
the inner rim of the disc, and the spiral structures at r>3a 
in the main body of the disc. The main torque frequencies 
must come from the interaction between the forcing binary 
quadrupolar potential and such overdensities propagating in 
the disc. To test this hypothesis, we mimic the situation by 
placing test particles in circular orbits at r — 2a and r = 
3.5a around a circular binary. We compute the torques and 




f[Po-] 

Figure 11. Simple test model for the periodicity generated in the 
outside disc, i.e. at r > a. In the top panel we show the temporal 
evolution of the torque, whereas in the bottom panel we show 
the Fourier transform. In each panel, the blue curve is obtained 
by placing three point masses at distance 1.6a, 2.1a, 3.5a from a 
circular binary, and the the red is the torque found in the isolO 
simulation. 

show them as blue lines in Fig. 1111 The natural frequency 
between a quadrupolar potential oscillating with frequency 
/1 and an overdensity orbiting at frequency / 2 is the beat 
frequency 2/i — f%, and this is what we see in the Fourier 
spectrum. The beats between the binary and the particles 
at r — 2a and at r — 3.5a give rise to sharp peaks seen 
at / « 1.35P _1 and / ~ 1.7P _1 , respectively. The fact 
that such peaks are much sharper in the iso simulations 
stems from two facts: (i) the binary does not significantly 
shrink in the iso runs, limiting the broadening of the spectral 
feature and (ii) the disc features (disc edge and spiral arms) 
are much more defined in this case, and the torques are 
well localized. This can be also seen by comparing the much 
neater torque structure in the bottom panel of Fig. |S]to the 
blurry one of Fig. [9] The / ~ P^ 1 peak of the blue line 
in Fig. [TT] is obtained by placing a third particle at r = 
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1.6a, a separation corresponding to the maximum negative 
torque inside the cavity. However, Fig. [7| shows that the 
peak at / w Pq" 1 is stronger at r > 1.8a than in the a < 
r < 1.8a region, meaning that it must be mostly directly 
related to the periodicity at which the binary rips gas off 
the disc edge (i.e. the binary period), and not to a specific 
radius in the cavity. Finally, as already noticed, torques in 
the binary region are related to the mass fuelling the two 
BHs, and therefore show periodicities related to the binary 
(/ = P _1 ), to the inner rim of the disc (/ ~ 0.25P ~ ) and 
to the beat between the two (/ « 0.75P _1 ). Note that the 
latter two are much more evident in the adia runs, where 
the disc rim overdensities feeding the accretion streams are 
more pronounced. 



5 ACCRETION 

In Tab. [2] we also show the average mass accretion rates. 
We recall here that a particle is considered accreted when 
it crosses the sink radius of one of the BHs. Then both its 
mass and its linear momentum are added to the BH, en- 
suring the global conservation of both quantities. The mass 
accretion rate is almost independent on the sink radius in 
the iso simulation, where instreaming gas settles in well de- 
fined accretion discs progressively loosing angular momen- 
tum. Conversely, in the adia case, the accretion rate scales 
almost linearly with the sink radius, suggesting a direct re- 
lation to the BH cross section (defined by the sink radius 
itself). 

It is particularly interesting to compare the accretion 
rate to the mass flow rates. Firstly we notice that in all 
simulations the accretion rate and the mass flow rates at 
r = a and r — 1.5a are nearly the same, meaning that the 
system is in a steady state configuration. Interestingly, the 
mass accretion rate is ~ 25% lower than the inflow rate 
at r = a. This means that not all the material crossing 
r — a is accreted, an assumption commonly adopted in grid 
simulations where the r < a region in excised. Part of the 
gas suffers a slingshot impacting back to the disc, affecting 
further the disc-binary angular momentum transfer. 

In the top panel of Fig. [12] we show the temporal evo- 
lution of the accretion rate on the primary hole, Mi , on the 
secondar y hole, Mi, an d the s u m of the two, M . As already 
found in ICuadra et all l|2009l ): iRoedig et all l|201ll h Mi is 
a factor ~ 2 larger than Mi, due to the closer interaction 
between the secondary BH and the disc. The relative mass 
growth 8M/M is therefore about six time faster for the sec- 
ondary BH, implying the binary mass ratio will tend toward 
q m 1 if this condition persists over the entire evolution of 
the system. Accretion rates are strongly modulated in time, 
and notably, M2 shows a striking periodicity related to the 
binary orbital period, whereas Mi is dominated by longer 
timescale fluctuation, related to the overdensities developing 
at the inner rim of the circumbinary disc. 

The Fourier transform of M is given in the lower panel 
of Fig. 1121 We observe the usual three distinctive peaks ob- 
served in the torques coming from the binary region, as 
shown by the central panel of the figure. It is clear that the 
accretion and inner-torque periodicities are intimately con- 
nected, both reflecting the temporal evolution of the mass 
inflow in the binary region (ii). 




Figure 12. Accretion onto the two BHs. Top box: accretion rate 
as a function of time. The red line is the total accretion rate while 
the blue and black lines are the accretion rates onto the secondary 
and the primary respectively. Lower box: Fourier transform of 
the total accretion rate (lower panel) and of the torques (upper 
panel). In the upper panel, the green line is the total torque (high- 
est peak normalized to 1) and the red line is the torque exerted 
by the material at r < a (multiplied by ten). 



Finally we can check the orientation of the an- 
gular momentum vector of the accreted material in 
the reference frame of the accreting BH. This num- 
ber quantifies the level of 'coherence' of the accretion 
flow, and has important consequences for the individ- 
ual BH spin evolution and the magnitude of gravita- 
tional recoil at the coalescence (Bogdanovic et all 2007; 



2010; 



Dotti et alj|20ld: iKesden et alj|20ld: IVolonteri et al.l 
Lousto fc Zlochowerl l201ll ; lLousto et al.l 120121 )' At each 
snapshot we therefore compute the average L z and L of 
the accreted particles in the accreting BH reference frame, 
and the angle 9 — arccos (L 2 /L) defining the degree of mis- 
alignment with respect to t he z axis. O u r resu lts are similar 
to what already found by iDotti et al.l (|2009l ). although in 
the cited study the spin evolution was studied only before 
the gas scouring. More specifically, an isothermal EoS leads 
to extremely coherent accretion flows, with fluctuations in 
the angular momentum direction of the mini-discs of few de- 
grees. In the adiabatic case, the orientation of the minidisc 
orbiting the secondary changes by at most « 20 degrees, and 
has an average excursion of 7 degrees, 50% bigger than the 
oscillations in the primary disc. Assuming that the BH spins 
have efficiently aligne d with the circum binary disc before the 
opening of the cavity (jPotti et al.ll201ol ) this implies that the 
efficient spin-up of the binary will continue following cavity 
opening. 



12 C. Roedig et al. 



6 COMPARISON WITH PREVIOUS WORKS 

The binary evolution and torque structure have been pre- 
viously investigated by MM08 and Sll, as described in the 
introduction. In particular, we can compare results regard- 
ing the average gravitational torques and the mass accretion 
rates. While both of them can be expressed in several ways, 
we find it practical to normalize these quantities to the disc 
mass, Md, initial binary period, Po, and initial binary angu- 
lar momentum magnitude, Lq. In these units we can write: 



M = IWdPo" 1 



(8) 
(9) 



Here Md = TrEpJ-p is the disc mass computed using the peak 
surface densitjjj E p . In both MM08 and Sll, r p ss 3a, close 
to what we find in our discs (see Fig. I10[) . Before proceed- 
ing with this comparison we should mention several issues 
that has to be considered. Firstly L = (M /A)(GM a ) 1/2 
for the circular equal mass binary simulated by MM08 and 
Sll, whereas in our simulation with q — 1/3 we have 
Lo = (3Mo/16)(GMoao) 1/ ^ 2 (assuming a circular binary). 
Secondly, the concept of 'accretion rate' is somewhat ill de- 
fined in grid based simulation with an excised central region. 
The inner boundary of the calculation is at r — a for MM08 
and r = 0.8a for Sll, and in both cases, whatever crosses 
that boundary is considered accreted. As already discussed 
(see Tab. [2]) ~ 25% of the mass crossing r = a is acceler- 
ated back to the disc edge extracting angular momentum 
from the binary (an effect already seen and described exten- 
sively by MM08 and Sll for the material inflowing in the 
cavity but not reaching the excised region). We can there- 
fore consider MM08 and Sll to overestimate the accretion 
rate by a factor of 25%. Lastly, MM08 and Sll can com- 
pute gravitational torques only for r > a, outside the binary 
corotation radiufl As shown by both studies (and indepen- 
dently recovered by us), gravitational torques exerted by the 
disc elements on the binary are negative in this region. 

For the sake of the comparison we use the iso05 and the 
adia05 runs. In the torque computation, we find Ti = —3.5 x 
1CT 3 for the former and Fi = -5.3 x 10~ 3 for the latter. This 
can be compared to -1.26 x 10" 3 of MM08 and -1.8 x 1CT 2 
of Sll. Our gravity torques are a factor 3-4 larger than 
MM08 and a factor 4-5 smaller than Sll. Concerning the 
accretion, we get Y2 = 1.1 x 1CP 3 for the iso05 run and 
T 2 = 6 x 10" 4 for the adia05 run, to be compared to Y2 = 
1.2 x 10" 4 of MM08 and T 2 = 6x 1(T 3 for Sll. Our accretion 
rates are a factor 5-10 larger than MM08 and a factor 5-10 
smaller than Sll. The similar scaling between (Tm i00 ]) and 
M is not surprising, since, as discussed in Section [4] the 
former is directly linked to the streams fuelling the BHs. 



7 DISCUSSION AND CONCLUSION 

The evolution of BHBs in circumbinary discs remains largely 
an open issue. Here we carried out a detailed study of 



8 Note that we have a physical disc mass in our simulations, but 
we use this quantity to ease the comparison with MM08 and Sll. 

9 Sll already pointed out the change of sign of the torques at 



the torques acting on the BHB, including the effect of 
the outer disc, the mass flowing in the cavity and, for 
the first time, of the gas entering the BHB region and 
eventually accreting onto the two BHs. We paid partic- 
ular attention to investigating the individual contribu- 
tion of different disc regions and we aimed at separat- 
ing the effect of accretion from the effect of gravitational 
torques. The emerging picture is far more complex than 
the often adopted simple assumption that resonant torques 
arising at the disc inner edge transfer angular momen- 
tum fro m the BHB to the disc, causing the orbital de- 
cay ( e. g. iPapaloizou fc Pringleill977l;lGoldreich fe Tremaind 
11980; 



Ivanov et al. 19991 ; Armitage fc Nataraian 20021 ; 



lHaiman et al. 2009) 



Firstly, as already pointed out by a number of authors 
(e.g. MM08, Sll), the presence of a BHB clearing a cavity in 
the disc does not prevent gas inflows and eventual accretion 
onto the two BHs. Such inflowing gas has a non negligible 
role in defining the BHB-disc mutual evolution; it forms 
streams that follow the BHs, eventually exerting a net neg- 
ative torque in the cavity region. The inflows catch-up with 
the BHs at a super-Keplerian speed, bend in front of them, 
thus exerting a net positive torque inside the binary corota- 
tion radius (i.e. at r < a). We therefore conclude that the 
origin of the dominating gravitational torques on the binary 
is purely kinematic, and can be understood without directly 
invoking resonances of the forcing binary potential with the 
disc. 

The strong positive gravitational torque is related to 
the gas inflows only, and not to the formation of minidiscs 
around the two BHs nor to the eventual accretion. This be- 
comes clear by inspecting Figs. [6] [10] and Tab. [2] although 
there is a factor up to ten difference in the mass bound to 
the BHs in minidiscs, and a factor of three difference in the 
accretion rate among the simulations, the positive torque is 
almost the same. This is because the gravitational torque is 
only related to the inflowing mass bending in front of the 
BHs, which is of the same order in all simulations. Note that 
different thermodynamics lead to different negative/positive 
torque balance. Although this might be an artefact of the 
sudden change of EoS in the iso simulations, such a balance 
has to be better understood in order to assess the fate of the 
binary. 

Torques exerted by gas in the main body of the disc 
(r > 2a) are instead negligible for the long term evolution 
of the system, since they average to zero (Fig. [p]F°l. Local 
torque maxima can be assoc iated to the 3:1 and 4:1 OLRs 
|Artvmowicz fc Lubow|[l994h . Given the nature of the forc- 
ing potential, a quadrupolar torque pattern naturally devel- 
ops (Figs. [H] and [§]) causing a time oscillating torque (Fig.Q) 
with characteristic frequencies related to the beat between 
twice the binary frequency and the characteristic orbital fre- 
quencies of inhomogeneities developing in the disc in form 
of inner edge overdensities and spiral arms (Fig. Ill[) . 

Accretion plays an important role in the binary angu- 
lar momentum budget. The binary can actually gain angular 
momentum even if it shrinks and its eccentricity increases. 
This is because accretion deposits a significant amount of an- 
gular momentum in the system by increasing its mass and 



Note that this region excludes the disc edge 
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mass ratio (see Eq. ([5])). The peculiar case of iso05 is quite 
interesting: here the total gravitational torque exerted by the 
gas is almost null (Fig. [6]). However, the BHB-disc mutual 
torques are responsible for the eccentricity growth, mean- 
ing that, if angular momentum was transferred by gravita- 
tional torques only, the BHB would be forced to expand. 
This is what we see in the upper-right panel of Fig. [3l where 
we show that the accretion-driven shrinking is faster than 
the actual evolution found in the simulation. In this case, 
gravitational torques alone would force the BHB to expand, 
whereas accretion is responsible for the shrinkage. Although 
whether this happens in Nature is questionable, it highlights 
the complexity of the BHB-disc interplay. 

We should, nonetheless, be careful in drawing any con- 
clusion about accretion from ou r simulations. If we scale our 
system to the fiducial binary of lRoedig et al.l (|201ll ) and set 
Mi = 2.6 x 10 6 M and a = 0.039 pc, then the accre- 
tion rate we find corresponds to about 20-40MEdd for the 
secondary BH, and 4-7MEdd for the primary BH. In such sit- 
uation, it is likely that most of the gas which is numerically 
accreted in our simulations will be expelled through outflows 
and winds. In this case, if the gas binds to the BHs before 
being expelled, it transfers its linear (an d angular) momen - 
tum to the holes, even if not accreted (|Nixon et al.ll20l"lT ). 
However, its mass does not add-up to the binary, making 
the question of angular momentum transfer more delicate 
and worthy of deeper and more refined investigations. 

The relatively low computational cost of our simula- 
tions allowed us to investigate different EoSs and sink radii 
''sink- The numerical size of r s i n k has a minor impact on the 
general behaviour of the system, affecting the mass bound 
to the BH in minidiscs and the accretion rate in the adia 
case. Conversely, the adopted EoS has a major impact on 
the results. In the iso runs, tenuous streams leak from an 
almost circular disc edge and fuel two well defined mini- 
discs, whereas streaming activity is more violent in the adia 
case, with thick spirals leaking from a deeply distorted disc 
edge, forming a diluted, oo-shaped cloud around the binary. 
This difference in streaming activity affects both the overall 
structure of the outer disc and the long-term binary evo- 
lution (as mentioned above). Although the extrapolation of 
the results of all our four simulations would imply a fast 
BHB coalescence (in less than 10 7 yr, for the fiducial binary 
parameters considered above), the almost perfect cancella- 
tion of the net gravitational torque in the iso simulations 
suggest that more realistic models for gas thermodynamics 
will be necessary to make clear-cut statements on this issue. 
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